Instability-induced fermion production in quantum field theory 
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Nonequilibrium instabilities are known to lead to exponential amplification of boson occupation 
numbers for low momentum modes on time scales much shorter than the asymptotic thermal equili- 
bration time. We show for Yukawa-type interactions that this growth induces very efficient fermion 
production, which proceeds with the maximum primary boson growth rate. The description is based 
on a 1/N expansion of the 2PI effective action to NLO including boson- fermion loops, which are 
crucial to observe this phenomenon. For long enough amplification in the boson sector, fermion 
production terminates when the thermal occupancy is reached in the infrared. At higher momenta, 
where boson occupation numbers are low, the fermion modes exhibit a power-law regime with ex- 
ponent two. 



INTRODUCTION 



Fermion production in the presence of large boson oc- 
cupation numbers is an important question for various 
aspects of early-universe cosmology as well as collision 
experiments of heavy nuclei in the laboratory. Large bo- 
son occupation numbers can lead to substantial devia- 
tions from perturbativc production processes. An impor- 
tant application concerns the question of thermalization 
following nonequilibrium instabilities in the context of 
early- universe infiaton reheating dynamics [l], H, H, 0, H[ , 
or in quantum chromodynamics (QCD) relevant for rel- 
ativistic heavy-ion collisions [1,0,1, IM E3- Nonequi- 
librium instabilities lead to exponential growth of boson 
occupation numbers in long wavelength modes on time 
scales much shorter than the asymptotic thermal equili- 
bration time. This is followed by a period during which 
bosons develop power cascades in a manner reminiscent 
of Kolmogorov wave turbulence [H], [HI |T3, [[H, . For 
long wavelengths a new class of nonperturbative scal- 
ing solutions has recently been found in the context of 
early- universe reheating dynamics [l3| . The correspond- 
ing nonthermal infrared fixed points [l3l , [l7j have the 
dramatic consequence that a diverging time scale exists 
far from equilibrium, which can prevent or substantially 
delay thermalization. Understanding these phenomena 
in the presence of fermions is an important step towards 
more realistic phenomenological applications. 

In this work we study nonequilibrium fermion pro- 
duction in Yukawa-type models following a spin- 
odal/tachyonic instability. More precisely, we consider 
a SU(2)l x SU(2)r symmetric theory for massless Dirac 
fermions coupled to a N = 4 vector of scalars in 3 + 1 
dimensions. The model incorporates the chiral symme- 
try of massless two-flavor QCD. It can also be consid- 
ered as a model for scalar infiaton dynamics coupled to 
fermions. For our purposes, a major aspect of this model 
is that controlled approximation schemes are available, 
which can describe the far-from-cquilibrium dynamics at 
early times as well as the late-time approach to thermal 
equilibrium. It is based on a 1 /N expansion to next-to- 



lcading order (NLO) of the two-particle irreducible (2PI) 
effective action out of equilibrium [l8|, [l9[ . The fermion- 
boson loop-contributions included in this approximation 
allow us for the first time to describe the phenomenon 
of instability-induced fermion production [49( : This pro- 
duction mechanism is based on the exponential growth 
of long-wavelength boson occupation numbers following 
a nonequilibrium instability. It is very efficient in the 
sense that the amplification of fermion occupation num- 
bers proceeds with the maximum primary growth rate 
of boson occupancies. This can lead to a fast approach 
to thermal equilibrium fermion distributions for low mo- 
mentum modes, though the bosons are still far from equi- 
librium at this stage. At higher momenta, where boson 
occupation numbers are low, the fermion occupancies ex- 
hibit a power-law regime with exponent two. Apart from 
a numerical solution without further assumptions, we 
present approximate analytical solutions which provide 
a detailed understanding of the underlying dynamics. 

Far-from-equilibrium fermion production is typically 
studied by integrating the Dirac equation in the pres- 
ence of a classical macroscopic field, such as in Refs. [22], 
H, H, [H, [13, H, M, [H| for scalar infiaton dy- 
namics coupled to fermions, or as in Ref. [12] for the 
production of quark pairs from classical gluon fields in 
the context of QCD. Fermion dynamics in the presence 
of inhomogcncous classical fields has been studied in an 
abelian Higgs model in 1 + 1 dimensions [33[. Analyti- 
cal studies of perturbativc dissipation and noise kernels 
for Yukawa- type models are presented in Ref. [11] . Stud- 
ies including quantum fluctuations, such that thermaliza- 
tion with the approach to Bose-Einstein and Fermi-Dirac 
distributions can be observed, have been performed in 
Refs. [l!| [H, [36[ for the model considered in this work. 
Here we extend these studies by including the NLO cor- 
rections in the 2PI 1/N expansion and apply the ap- 
proximation to the physics of nonequilibrium instabili- 
ties. This has not been done before including fermion 
fluctuations. In contrast, the bosonic part of our model 
has been studied in great detail and the 2PI 1/N expan- 
sion to NLO is known to provide a quantitative descrip- 
tion of nonpcrturbatively large boson fluctuations for the 



2 



considered case with TV = 4 [f| [HI, Here we eval- 
uate the 2PI effective action for vanishing macroscopic 
field. The observed phenomena in the fermionic sector 
are, therefore, entirely due to quantum corrections that 
are neglected in previous studies where the Dirac equa- 
tion is solved in the presence of a classical macroscopic 
field. In particular, our findings of an induced exponen- 
tial growth of fermion occupation number is distinct from 
the mechanism of resonant fermion production due to 
an oscillating inflaton field |24[. Similarly, the observed 
power-law behavior at high momenta is shown to be a 
consequence of quantum fluctuations beyond mass-like 
corrections, which are absent for the considered massless 
fermions because of chiral symmetry. It will be an im- 
portant further step beyond our current work to include 
non- vanishing macroscopic field values, in order to esti- 
mate the importance of quantum corrections to the Dirac 
equation in the presence of classical fields. 

The paper is organized as follows. In the next sec- 
tion we introduce our model with its approximation and 
initial conditions. In Sec. Illll thc boson dynamics is dis- 
cussed, where we will describe the early-time behavior 
following the noncquilibrium instability as well as the 
subsequent evolution to a nonthermal fixed point. Sec. lIVI 
is dedicated to a detailed analysis of the fermion dynam- 
ics, with the induced fermion production and the emer- 
gence of a power-law regime. In Sec. |V] we summarize 
and present our conclusions. There are three Appendices. 
In the first we give numerical and analytical results for 
the time evolution of the spectral function. Appendix [B] 
presents some details about the chiral decomposition of 
the fermion two-point correlation functions, and the last 
Appendix is about the employed numerical method. 



a metric with components r)ui> — diag(l, — 1, — 1, — 1), 
where summation over repeated indices is implied. The 
fermion fields i/j(x) with x = {a;°,x} and %p = 1/^7° 
are coupled via the Yukawa interaction g to the scalars 
d> = {cr, 7T 1 , 7r 2 , 7r 3 }, whose self-coupling is A and mass 
parameter m. 

Far-from-equilibrium dynamics as well as subsequent 
thermalization can be efficiently studied in this model 
using a resummed 1/N expansion beyond leading order 
(l8j |. based on the two-particle irreducible (2PI) effective 
action |1, S EE EH El ■ A major benefit of the 1/N 
approximation is the possibility to quantitatively address 
dynamics even in those cases where fluctuations become 
nonperturbative in the coupling, which will be relevant 
for the physics considered in this work. In noncquilibrium 
quantum field theory the 2PI resummed 1/iV expansion 
is required to cure the secularity problem of the stan- 
dard (IPI) 1/N expansion, which allows one to study 
the late-time behavior of quantum fields [Hj]. For the 
quantum theory with action (JTJ) the relevant nonequilib- 
rium evolution equations have been derived in Ref. (l9| 
to which we refer for further details. Here we summarize 
the derivation and give the equations which we use for 
our analysis. 

The most general 2PI effective action can be written as 
a functional of one-point and two-point correlation func- 
tions, such as (4>(x)) and (<f>(x)4>(y)) for scalars and sim- 
ilarly with fermion fields. Here the brackets (. . .) denote 
traces over an initial density matrix, which is discussed 
below. We will always evaluate the functional for van- 
ishing one-point functions, i.e. (0) = (ifi) = (i/)) = 0, and 
denote the nonvanishing two-point correlation functions 
by G(x,y) = {<t>(x)<t>(y)) and D(x,y) = {i>(x)$Ml In 
this case the 2PI effective action can be written as I38L H 



II. MODEL AND APPROXIMATION 

A. SU(2) L x SU(2) R symmetric model 

We consider a class of quantum field theoretical mod- 
els for fermions interacting with scalar bosons. This 
type of models has been widely used in very different 
contexts, such as fermionic preheating after inflation or 
as effective theories for low-energy hadron physics in 
quantum chromodynamics (QCD). Here we consider a 
SU(2) L x SU{2) R ~ 0(4) symmetric theory with JV> = 2 
species of massless Dirac fermions coupled to N s = 4 
scalars. The classical action reads 



S = I d 4 x ^ ipij^d^ip + j^r- ip (c + 1757-' 7r') ip 



-m 2 (cr 2 + 7rV) 

2 V ; 4!A S 



(o» + ,(1) 



where r (I = 1, 2, 3) are the standard Pauli matrices and 
7 M (/i = 0, 1, 2, 3) denote the Dirac matrices. We employ 



T[G,D] = -tTr(D^D) + -Tr(G^G) 



iTrln(Z)- 1 ) + -Trln (G -1 ) 



r 2 [G,L>] + const. (2) 

with the free scalar inverse propagator 

*G^ ab {x, y) = - (d^ + m 2 ) 5^ (x - y) S ab (3) 

for a, b = 1, . . . , N s and the massless free inverse fermion 
propagator 



iD-^{x,y) - i-fd^\x-y)8i. 



(4) 



for i,j 



1. 



,N f 



The corresponding components 



of the full boson (fermion) propagator are denoted by 
G a b(x,y) (Dij(x,y)). The traces in indicate inte- 
gration over time along the closed time path [43| and 
over three spatial coordinates as well as summation over 
internal indices. The functional T2[G, D] only contains 
contributions from two-particle irreducible diagrams with 
propagator lines associated to G a b{x,y) and Dij(x,y) 



The equations of motion for the full propagators are 
obtained from the stationarity conditions [3a . [39( 



ST[G,D] 
5Dij(x,y) 



0. 



5T[G, D] 
SG ab (x,y) 



0. 



(5) 



where we omit Dirac indices in the notation. Self-energies 
arc defined by functional differentiation of the 2PI part 
T 2 [D,G] as 



2i 



: 6T 2 [G, D] 
' 8Dji{y,x) 
; ST 2 [G,D] 
8G ab (x 7 y) ' 



(6) 



Because of SU(2)l x SU(2)r ~ 0(4) symmetry one can 
always choose the propagators to be diagonal, such that 
G a b = G5 a b, Dij = D Sij and similarly for the self- 
energies. 

It is convenient to decompose the time-ordered two- 
point correlation functions D(x 1 y) and G(x, y) into a real 
and an imaginary part. For the fermion correlator we 
use [H] 

D(x, y) = F(x, y) - % -p{x, y) sign(a-° - y°) , (8) 

with the real statistical two-point function F(x, y) and 
the spectral function p(x,y). For fermions F is deter- 
mined by the expectation value of the commutator of 
two fermion field operators, while p is described in terms 
of the anti-commutator of two fields [42j ■ A similar de- 
composition can be done for the boson correlator with 

G(x, y) = F^(x, y) - -p^x, y) sign(a; - y°) , (9) 

where, in contrast to the fermion case, F^ is determined 
by the expectation value of the anti-commutator of two 
boson field operators and the spectral function p^ by 
their commutator [42| . Correspondingly, the decompo- 
sition identities for the self-energies © and ([7]) read 

S(x,y) = C(x,y)- l -A{x,y)sign{x° - y°) , (10) 
11(1, ») = -iIl local (x)6^(x-y) 

+ n F (x, y) - l -U p (x, y) sign^ - y°) . (11) 

Here we take into account a possible space-time depen- 
dent local term IIi oca i(a;) for the boson self-energy, which 
is not present for the fermion self-energy [lj|. 

We will consider spatially homogeneous and isotropic 
systems. The chiral symmetry of the model ([1) together 
with parity and CP invariancc imply that for the fermion 
two-point functions only the vector components 

K = |tr( 7 ^F) , p v = Jtr( 7 V) , (12) 
C v = \t r (rC) , A^r = Jtr( 7 M) (13) 



are non- vanishing [19J, where the trace acts in Dirac 
space. The vector components can be decomposed in 
spatial Fourier space as 



F 1 ' 



(x,y) = f e^^F^x ,y°;p), (14) 
J p 



where we use the abbreviation J p = J d 3 p/(2ir) 3 , with 



F°(x ,y° ]P ) = F v (x°,y°;\p\), 
F v (x°,y°;p) = £-F v (x°,y°;\p 



(15) 



and equivalently for the other two-point functions. 

Denoting times as x° = t, y° = t' etc. and setting 
the initial time to zero (5lj . the equations of motion ([5]) 
for the fermion two-point functions F v , Fy, p v and pv 
read [H: 

id t F v (t,t';\p\) = \p\F v (t,t';\p\) 

t 

+ J dt"{A v (t,t";\p\)F v (t",t';\p\) 
o 

-A v (t,t";\p\)F v (t",t';\p\)} 
i' 

dt"{C v (t,t"-Ap\)Pv(t",t';\p\) 

o 

-C v (t,t";\p\)p v (t",t';\p\)}, (16) 

id t F v (t,t';\p\) = \p\F v (t,t';\p\) 
i 

+ J dt"{A v (t,t";\p\)F v (t",t';\p\) 
o 

-A v (t,t";\p\)F v (t",t';\p\)} 



dt"{C v (t,t";\p\)p v (t",t';\p\) 



-C v {t,t"-\p\)p v (t",t'-\p\)} 



(17) 



id tPv {t,t';\p\) = |p|MM';IpI) 
+ J dt"{A v (tX;\p\) P °v(t",t'-Ap\) 

V 

-Av(t,t";\p\)p v (t",t';\p\)}, (18) 

id tPv (t,t';\p\) = \p\p° v (t,t';\p\) 
t 

+ J dt"{A v (t,t";\p\)p v (t",t';\p\) 
v 

-A v (tX;\p\)p v (t",t';\p\)} ■ (is) 

Similarly, the equations of motion for the boson two-point 



functions F^ and p^ are 

{d 2 + p 2 + Af 2 (i)}^(i,i';|p|) = 
/ 

-J dt"u p {t,t"-M)m"J-M) 

/' 

+ 1 dt"U F (t,t";\p\)p 4> (t",t';\p\), (20) 
o 

{a 2 + p 2 + A/ 2 (i)} P0 (i,t';|p|) = 
t 

-J dt"lL p (t,t"-APl)P^t",t';\p\), (21) 
/' 

with the mass-like term 

M 2 (t) = m 2 + II local (t), (22) 

which depends on time only because of the assumed spa- 
tial homogeneity. The above equations have to be regu- 
larized, which we will do by implementing a momentum 
cutoff following Ref. as described in Appendix [U] 

One observes the following symmetry or anti- 
symmetry relations with respect to interchanging time 
arguments [l9l ]: 

F v (t,t';\p\) = F v (t',t;\p\), 

F°(t,t';\p\) = -F°(t',t;\p\), 

p v (t,t';\p\) = -p v (t',t;\p\), 

P a v (t,t';\p\) = p° v (t',t;\p\). (23) 

For the Fourier coefficients of the boson two-point func- 
tions the corresponding relations are 

Ft(t,t';\p\) = F^t',t;\p\), 

P<i> (t,t';\p\) = -p</>(t',t;\p\). (24) 

B. Approximation 

For given T 2 [G, D] the self-energies are expressed in 
terms of self-consistently dressed propagators G and D 
according to ([6]) and ([7]). Here we consider the non- 
perturbative 2P1 1/N expansion to next-to- leading or- 
der (NLO) in the number of boson fields N s [l8]. For 
bosons the physics of noncquilibrium instabilities can- 
not be described by an expansion in the coupling A even 
for weak couplings, since occupation numbers can grow 
parametrically large of order 1/A @. In contrast, oc- 
cupation numbers for fermions are bounded due to the 
Pauli principle and a 2PI coupling expansion in powers of 
the Yukawa coupling g will be employed to describe the 
fermion fluctuations. This requires to sum an infinite se- 
ries of boson diagrams, which is displayed in Fig. [1] along 
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r r 



FIG. 1: Diagrammatic representation of the employed ap- 
proximation for T2. Solid lines are associated to self- 
consistently dressed boson propagators G and dashed lines to 
fermion propagators D. The dots indicate an infinite series 
of diagrams, each obtained from the previous one by adding 
another loop in the bubble ring. 



with the single boson-fermion diagram contributing at 
this order. 

The corresponding self-energies entering the evolution 
equations (QI|)-(|2l]) are for N s = 4 and N f = 2 

C£(t,i';p) = -g 2 ^{F(^,i';q)^(M';p-q) 

-^(t,l/;q)p^t,t';p-q)\ , (25) 

A£(*,t';p) = -g 2 y" j^(M';q)p0(M';p-q) 

4y;(i,t';q)F (M';p-q)} (26) 
for the fermion part and 

U F (t,t';p) =~^J |i^(M';p-q) 
x I F (t, t'; q) - ip (t, t'; p - q)J„(t, t'; q)| 
-2. 9 2 J ^F$(t, t'; q)F V)M (i, t'; p q) 
~p^{t,t'^)p v ^(t,t'-,p-q)X, (27) 

n p (M';p) = ~ J ^(t.t'sp-q) 

x I p (t, t'; q) + p^t, t'-p- q)I F (t, t'; q) J 

-V J \p^{t, t'; q)F Vtfl (t, t'; p - q)| (28) 

for the boson part. Here the functions I F and I p con- 
tain the summation of the infinite series of boson loop 
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diagrams at NLO: 



t 

~P4,(t, t'; p - q>*(t, t'; q) - J dt"l p (t, t"; p) 

o 

x (f^", f- P - q)W, q) - jM*". p - q) 

xp^(i",t';q)J +2 y" dt'7 F (t,t";p)F (t",i';p-q) 



xp*(t",f;q) 



(29) 



I P (t,t';p) = ^ j |i^(t,*';p-q)p^(t,t';q) 

- ^ df7 / ,(t,*";p)^(t ,/ ,t';p-q)^(f ,/ ,t , ;q)|. (30) 



Finally, the mass-like term (|22j) reads 



A 



M^) = m^ + - / F (t,i;q) 



(31) 



This approximation leads to a closed set of equations for 
the two-point functions Fy, Fy, p v , py, F,p and p^. 



C. Initial conditions 

In order to solve the (integro-) differential equations 
of motion (fT6]) - (f2~Tj) with the approximation (|2"5|) - (f3Tj) . 
we have to specify initial conditions. The initial condi- 
tions for the spectral functions are fixed by their (anti-) 
commutation relations i.e. they are given by the 

equal-time properties for the fermion spectral function 

p° v (t,t;p) = i, p v (t,t;p) = (32) 

and for the boson spectral function 

P<t>{t,t\p) = 0, d t P4,{t,t' ;p)\ t=v = 1, (33) 

which are valid for all times t. For the statistical func- 
tions we parametrize the Gaussian initial conditions as 

F v (t,t';p)\ t=t , =0 = 0, (34) 
F v (t,t';p)\ t=t>=Q = 1 (35) 



for the fermions and 



F4,{t,t';p)\ t =v=t 



1 



(36) 



2eo(p) ' 

d t F^t,t';p)\ t=t , =0 = 0, (37) 
£ o(p) 



dtdfF^tJ;?) 



t=t'=0 



(38) 



for the bosons. Here eo(p) = e(t = 0;p) plays the role 
of a quasi-particle energy at initial time, which we define 
using [TU 



e(t,p) = 



F^(t, t';p) 



(39) 



We note that the equal-time correlation function 
F®(t, t;\p\) corresponds to the conserved net charge den- 
sity [19| . which vanishes identically at all times for our 
initial conditions. We will consider the time-dependent 
quantities 



n v ,(t, p) 
n^>(t, p) 



1 



-F v (t,t;p), 



F^{t,t;p)e{t,p) 



(40) 
(41) 



which may be associated to effective fermion and boson 
occupation numbers, respectively. 

In our simulations we start from a vacuum-like state 
with n^{t = 0, p) = and n^it = 0, p) = 0. We induce a 
spinodal (" tachyonic" ) instability with a so-called quench 
by a sudden change in the sign in front of the m 2 term 
appearing in (|31| . More precisely, we employ 



(42) 



M 2 (t = 0) = m 2 , 



M 2 (t >0) = -m 2 + - I ^F (t,t;q) 



where the constant subtraction of quadratically divergent 
terms in the integrand of (|4"2")) ensures M 2 (t = + ) = 
— m 2 for the employed vacuum-like initial condition (|36|) . 
We choose m as our scale and express all quantities in 
appropriate powers of it. 

Though all results that will be shown are obtained from 
the above initial conditions, we note that the evolution 
quickly becomes rather insensitive to the details at initial 
time. We will see that the spinodal instability leads to 
an exponential growth of fluctuations which soon domi- 
nate over the initial values of correlators. Late stages in 
the far-from-equilibrium evolution can even become uni- 
versal in the sense that characteristic properties become 
insensitive to the details of the underlying microscopic 
field theory, which will be explained in the following. 



III. BOSON DYNAMICS 
A. Spinodal/tachyonic instability 

Since the evolution starts from a vacuum-like state, at 
very early times almost no fluctuations are present for 
weak couplings. Accordingly, neglecting the self-energy 
corrections in (|20[) and (f2"2"]) the initial evolution for the 
Fourier modes of the boson statistical propagator may be 
approximately described by 



[d 2 + p 2 -m 2 ]F (M';|p|) ~ 0. 



(43) 
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FIG. 2: Occupation numbers of boson modes as a function 
of time for given boson self-coupling A and different values of 
the Yukawa coupling g to fermions. The upper graph shows 
the results for the purely bosonic theory with no fermions. 
g = (dashed lines), compared to a weakly coupled theory 
with g\ // 2 = 0.1 (solid lines). The lower graph shows the 
behavior for gy/2 = 1. Here and in all following figures the 
quantities are given in units of the mass parameter m. 

Because of the negative sign in front of the m 2 > term, 
the modes with |p| < m exhibit exponential growth and 
the dominant solution behaves as 

*V(m';p) = A)e 7(p) {t+t '\ (44) 

where the amplitude Aq has the dimension of an inverse 
mass and we approximately have A^m ~ 1/4 for our ini- 
tial conditions. The momentum dependent growth rate 
is given by 

7(p) = yjm 2 - p 2 . (45) 

This classical instability leads to a strongest amplification 
with rate 70 = 7(p = 0) for the mode with the smallest 
momentum. This "primary" growth stage is character- 
ized by a limited range of momenta |p| < m for which 
amplification can be observed. During the subsequent 
exponential growth the fluctuations become larger and 




t 



FIG. 3: Time evolution of the bosonic statistical two-point 
function F r f,(x,x) for different self-couplings and g\pl = 10 -2 . 

the self-energy corrections in pp|) , (f!H]) and (J^U) become 
relevant. These self-energy corrections incorporate non- 
linear effects, such as scattering and off-shell dynamics. 
For the bosonic model, in the absence of fermions, this 
has been discussed in detail for parametric resonance as 
well as spinodal instabilities in Refs. 0,113, Ell- 
in the upper graph of Fig. [2] we show our results 
for the evolution of different boson occupation numbcr 
modes (|4"Tj) as a function of time for a scalar self-coupling 
A = 0.01. The graph compares two runs, one without 
fermions [g = 0) and the other with a weak Yukawa cou- 
pling of g\[2 = 0.1. On the logarithmic plot one observes 
that the boson occupation numbers are practically not 
affected by the weak coupling to the fermions. The time 
evolution exhibits the characteristic stages explained in 
detail in the purely bosonic studies of Refs. 0, E3, Ell- 
Low- momentum modes with |p| < m show primary ex- 
ponential growth, whereas higher momentum modes be- 
come exponentially amplified at a later ("secondary") 
stage with faster growth rates due to nonlinearities that 
were built up by primary modes. After the fast period 
of exponential growth the dynamics slows down consid- 
erably. At this stage the nonlinearities are nonpertur- 
batively large, i.e. paramctrically F$ ~ 0(1/ X) and all 
processes become of order unity. This is illustrated for 
F c j > (x, x) = J F$(t, t; p) for different values of the boson 
self-coupling in Fig. [3] One observes that the ratio of this 
quantity taken at two coupling values, A and A', is well 
described by 

F<j,{x,x)\\' A 

The time at which this parametric dependence can be 
first observed we call t$. This time was estimated in 
Ref. [11] for the case of parametric resonance and we can 
adopt the method. More precisely, we consider the time 
when the local self-energy contribution (" tadpole" ) of the 
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FIG. 4: Evolution of a low-momentum Fourier mode 
F$(t,t; |p|) for different couplings A and fixed g. The vertical 
dotted lines indicate the analytical estimates for the times at 
which the exponential growth period terminates: t$ = 7.2/m, 
8.5/m, 9.7/m, 11.0/m for A = 10" 2 , 10" 3 , 10" 4 , 10~ 5 . 



mass-like term (f3"Tj) becomes one: 



A^o e 27ot 
32( 7( 7 1 7rt) 3 / 2 



(47) 



For the evaluation of the momentum integral we took 
(|44[) and applied a saddle-point approximation around 
|p| =0. The last equality of (J47J) yields 



2 7 o V A 



470 



In 



7T^ 



,2/3 -1/3 
, ^0 /0 



(48) 



Plugging in Agm ~ 1/4 and using 70 = m one finds 
from the solutions of (J38J) that ~ 7.2/m for A = 10~ 2 , 
i ~ 8.5/m for A = 1(T 3 , ~ 9.7/m for A = 1CT 4 and 
t<j, ~ 11.0/m for A = 10~ 5 . Fig. [4]shows a comparison of 
the analytical estimates with the numerical solutions for 
the nonequilibrium time evolution. 

Occupation numbers as a function of momentum for 
different times are shown in Fig. The left column of 
that figure displays boson distributions for different cou- 
plings as indicated in the figure. The right column shows 
the corresponding behavior for the fermions, which will 
be discussed below in Sec. IIVI The time evolution of 
the boson occupation numbers is characterized by the 
exponential primary growth of low-momentum modes, 
and a subsequent broadening of the distributions as a 
consequence of enhanced secondary growth rates. Even 
though the Yukawa coupling changes one order of mag- 
nitude from gy/2 = 0.1 to g\[2 = 0.01, the corresponding 
graphs with scalar self-coupling A = 0.01 do not show 
any significant differences. This confirms that weakly 
coupled fermions have no significant affect on the boson 
dynamics. 

If only the boson dynamics is considered with- 
out fermions, alternative nonperturbative approximation 



schemes are available that can be used to test the 2PI 
1/N expansion [H, [3?], EE EH- Since the exponential 
growth is induced by an instability of the free-field-type 
equation (|43[) at sufficiently early times quantum correc- 
tions may be neglected. Furthermore, occupation num- 
bers become large at later times because of the exponen- 
tial growth and one expects an accurate description of 
the quantum dynamics using the classical-statistical field 
theory approximation for highly occupied modes. The 
classical-statistical simulation result is obtained from re- 
peatedly solving the classical field equation of motion 
numerically and Monte Carlo sampling of the initial con- 
ditions. Fig. [6] shows a comparison of the purely boson 
dynamics, represented by F<p(t,t; |p|), from the quantum 
2PI 1/N expansion and the classical-statistical simula- 
tion for the initial conditions (|3l)|) - (f3"8")) and A = 0.1. The 
level of agreement between both calculations is remark- 
able. For low momenta, where occupation numbers are 
high, the classical-statistical simulation is expected to be 
accurate up to controlled statistical errors and the com- 
parison points out the ability of the NLO approxima- 
tion to describe the strongly nonlinear dynamics. On the 
other hand, from the comparison one also observes that 
for the considered small coupling the quantum correc- 
tions are small even for lower-occupied modes at higher 
momenta for not too late times. For the asymptotic ap- 
proach to quantum thermal equilibrium characterized by 
a Bosc-Einstcin distribution the classical-statistical ap- 
proximation becomes invalid, while the quantum evolu- 
tion using the 2PI 1/N expansion to NLO approaches a 
Bose-Einstein distribution [lj, H3, EE El] ■ 



B. Subsequent power-law behavior 

The approach to thermal equilibrium can be substan- 
tially delayed in the presence of a nonequilibrium insta- 
bility [T2I ]. This happens because after the exponential 
growth period the evolution is driven towards a nonther- 
mal infrared fixed point. This was pointed out in Ref. fl3j 
for the case of a parametric resonance instability and 
discussed in terms of nonthermal renormalization group 
fixed points in Ref. [l7| . In order to demonstrate that 
a similar behavior also occurs starting from a spinodal 
or tachyonic instability, we present in Fig. [7] classical- 
statistical simulation results for much later times than in 
Fig. [SJ Shown are snapshots of the occupation number 
as a function of momenta for different times, where the 
black dots are taken at the latest time t = 4900/m. A 
straight line on the double-logarithmic plot corresponds 
to a power-law behavior, 



n<t,(t, |p|) 



1 



(49) 



where the occupation number exponent k at low- 
momenta is well approximated by kj r — 4 in agreement 
with the analytical estimates in Refs. [IEII3 f° r the non- 
thermal fixed point as well as the numerical results start- 
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FIG. 5: Occupation numbers for bosons (left) and fermions (right) as a function of momentum at different times. The graphs 
from top to bottom correspond to different couplings as indicated in the figures. The dotted lines in the graphs for the fermion 
distributions display power-law behavior n^,(]p|) ~ |p| with the exponent k^, = 2 for comparison (see Sec. IIV|I . 
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FIG. 6: Comparison between spectra of the boson correla- 
tion function F$(t,t; |p|) in a 2PI (solid lines) and a classical- 
statistical (dotted lines) simulation at different times for 
A = 0.1. The spectra for times larger than t = 5/m are 
multiplied by a factor of 10 in order to make separate lines 
better visible. 
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FIG. 7: The boson occupation numbers n$(t, |p|) from 
classical-statistical simulations for A = 0.1. Shown are in gray 
the distributions starting at time t = 900/m (from right), con- 
tinuing with a time step of At = 1000/m to the black dots 
which correspond to t = 4900/m. 



ing from parametric resonance [13| . At higher momenta, 
where paramctrically n^(t, |p|) < V^i one observes 
the transition to a classical-statistical thermal distribu- 
tion for which the occupation number exponent becomes 
kuv — 1- Towards higher momenta than those shown in 
Fig. [7] the distribution still drops rapidly. The classical- 
statistical theory will finally occupy all high-momentum 
modes according to a power-law with ktjv — 1, which cor- 
responds to the Rayleigh- Jeans ultraviolet divergence. In 
contrast, quantum corrections in the 2PI 1 /N approxima- 
tion prevent a power-law distribution at high momenta 
and lead to finite results [HI H H3, ES EH • 

If the Yukawa coupling g is increased in the pres- 
ence of fermions, the dynamics of the bosons becomes 
substantially affected and their behavior can no longer 
be understood separately. This is demonstrated in the 
lower graph of Fig. rj] for gy/2 = 1. We emphasize that 
for g f» 1 it may not be quantitatively justified to ne- 
glect higher loop-orders including fcrmion propagators 
beyond the two-loop graph taken into account in T2 as 
shown in Fig. [T] However, the observed increase of the 
growth rates comes from a negative mass-squared con- 
tribution for the boson self-energies, which is induced by 
fermion fluctuations. The very same mechanism oper- 
ates already for small couplings with g <C 1 and just 
the quantitative corrections are not sizeable enough in 
this case to be clearly visible in the figures. The induced 
negative mass-squared contribution is a well-known phe- 
nomenon, which has been analyzed in great detail also 
in the context of nonperturbative renormalization group 
studies of this model in thermal equilibrium [47j . Even 
though the boson dynamics is not much influenced by 
a weak coupling to fermions and sizeable effects require 
stronger couplings, a similar observation can not be made 
for the fermion dynamics. It turns out that the latter is 



always substantially affected by the highly occupied bo- 
son modes, since even for weak Yukawa couplings they 
induce an exponential growth in the fermion correlation 
functions as is explained in the following. 



IV. FERMION DYNAMICS 

A. Instability-induced fermion production 

In Fig. [H] different modes of the fermion occupation 
number (|40[) are presented as a function of time for two 
different Yukawa couplings, g\[2 = 0.1 (upper graph) 
and g\/2 = 1 (lower graph), for A = 0.01. The figure 
shows the behavior of the fermions, which corresponds 
to the evolution of the bosons presented in Fig. O Af- 
ter a short initial period one observes an exponential 
growth of fcrmion modes. The fast dynamics then slows 
down rather abruptly and approaches a quasi-stationary 
regime. A most characteristic property of the fast pe- 
riod is that the different momentum modes get exponen- 
tially amplified with approximately the same growth rate, 
which can be seen from the slopes of the curves in Fig. [5] 
for a wide momentum range |p|/m = 0.1 - 9.3. This is 
in contrast to the amplification of boson modes, where 
the primary growth rates show a significant momentum 
dependence according to (|45|) and the secondary rates 
equal multiples of the maximum primary growth rate (see 
Fig. [2]). Similar to the behavior of the boson secondaries, 
one observes that higher-momentum fermion modes are 
amplified at later times. Though fermion production oc- 
curs practically over the whole spectrum, this delay leads 
to a decrease of the amplitude of the occupation number 
modes with increasing momenta. One also notes that in 
our simulations n^it, p) never exceeds 1/2, despite the 
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FIG. 8: Occupation numbers of fermion modes as a function 
of time with couplings as indicated in the graphs. The ob- 
served growth rate of modes is approximately independent of 
momenta. In the lower graph only two modes are labeled, the 
three other modes correspond to the same momenta as in the 
upper graph. 



fact that its definition allows occupancies in the range 
< n^p{t, p) < 1 in accordance with the Pauli principle. 
We note that n,Fu(t, p = 0) = 1/2 corresponds to the 
value of the low-momentum distribution in thermal equi- 
librium, which is discussed in more detail in Sec. lIV Bl In 
the lower graph of Fig. [8] we show the same evolution but 
with a stronger Yukawa coupling gV2 = 1. One observes 
a more efficient fermion production over the spectrum, 
i.e. the decrease of the amplitude with higher momentum 
is less pronounced than for weak couplings. The growth 
rate is larger and the system enters the quasi-stationary 
evolution regime earlier. Corresponding plots for occu- 
pation numbers as a function of momentum for different 
times are shown in Fig. [5] above. The right column of 
that figure displays the fermion distributions for differ- 
ent couplings as indicated in the figure. The emergence 
of the characteristic power-law behavior for higher mo- 
menta seen in Fig. O will be discussed below in Sec. lIVBl 
In the following we derive an approximate analytical 
solution for the fermion dynamics, which explains the 



mechanism of the observed exponential growth and ma- 
jor characteristics. For weak Yukawa coupling the evo- 
lution of the boson correlation functions can be estab- 
lished without taking into account the fermion dynamics 
as discussed above. This is our starting point for an 
approximate evaluation of the fermion dynamics. We 
consider the equation of motion p7[) and concentrate 
on the equal-time correlator, Fy (t, t; p) , since this de- 
termines the behavior of the effective particle number 
PP)) . Because of the symmetry relations (f2"3|) the term in 
(jTTJ) proportional to the momentum |p| vanishes. With 
the definition of the self-energies (f2"5| and (f2T))) there are 
eight terms, each coming with a product of two fermion 
and one boson two-point correlation function. The self- 
energy contributions can be visualized by a loop diagram 
as shown in Fig. [5J where internal lines represent the var- 
ious propagators as indicated in the diagram. Because 
of the exponential growth of boson occupation numbers, 
the value of the respective equal-time correlation func- 
tion F^{t,t;p) becomes quickly many orders of magni- 
tude larger than any other function. In particular, it be- 
comes much larger than the spectral function p^t, t'; p), 
which even vanishes at equal times. (Sec Appendix [Al for 
a detailed discussion of the unequal-time properties of 
the spectral function.) This implies that all terms which 
include F^{t,t';p) will dominate the self-energy correc- 
tions in (jTTJ). The dominance is best realized if the boson 
self-coupling is small since the amplitude of F$ becomes 
at the end of the exponential growth period paramctri- 
cally of order 1/A according to (J3S1). Hence our analytical 
estimates will require small A. Neglecting the self-energy 
contributions without F^(t, t'; p) leads to four remaining 
terms and (TTj) becomes for equal times 

i{d t F v (t,t';\p\)}\ t >=t - 



g 2 dt" 



o 

_p__q 

IpI |q| 



p° v (t,t";\c l \)F 4> (t,t";\p- (l \)F v (t",t;\p\) 



M*,* ,/ ;|q|)^(*»*";|p-q|)^v(*".*;|p|) 



+ g 2 dt" 



o 

_p__q 

IpI |q| 



F v (t,t";h\)F 4> (t,t";\p-q\) P v(t",t;\p\) 



F v (t,t";\q\)F^t,t";\p^ci\)p a v (t",t;\p\) 



(50) 



Since we are considering the evolution of the equal- 
time correlation function, we employ for the derivative 
on the left hand side of (l50l): 



d t F v {t,t;p) = {d t ,F v {t\t-p)}\ t , =t 
+ {d t ,F v (t,t' ] p)}\ t ,= t 
= 2{d t >F v {t' 1 t-p)}\e =t , (51) 

where we used in the second equation that Fy(t , t] p) is 
symmetric under time exchange according to 
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FIG. 9: Topology of the fermion self-energy diagram. The 
internal lines represent the various propagators indicated in 
the diagram. 



For not too early times the integral in ([50)1 is domi- 
nated by momenta q ~ p for small enough p. To observe 
this we consider the primary growth rate (|45[) for p <C m, 
which is approximately given by 



7(P) 



2m 



(52) 



With (|44[) and 70 = m the low-momentum modes behave 



a.s 



F^t, i; p - q) ~ A e 2 ^ e-^P-q) 2 * , (53) 



such that the integrand of ([5(7)) is dominated by q ~ p 
for sufficiently large t. This holds also in the presence 
of additional powers of momenta in the integrand, such 
as coming from the measure of the integral. The sum 
appearing in the integrand of ([SITj) consists of pairs of 
terms, which are the same if q = p. Accordingly, the 
equation of motion may be approximated by 

t 

id t F v (t,t;\p\) ~ -4g 2 Jdt"J^(t,t";\p-q\) 



p» v (t, t"; \p\)F v (t",t; |p|) - p v (t,t"; |p|) IP 

(54) 

During the time interval of exponential growth of 
F^(i, i'; |p|) the latest times dominate the time- or 
" memory- integral" in f|54[) . The characteristic time scale 
for these dominant contributions is given by the primary 
growth rate 70. Followin g si milar calculations for bo- 
son dynamics in Refs. [El |42|. we will employ a mem- 
ory expansion for an approximate analytical description 
of the fermion evolution. Integrals over the entire past, 
J dt' , are replaced by integrals with a finite time range, 

ft-c/ta ^ e W ^ see tnat during the growth period 
this memory restriction leads to a good description of 
the 1 full numerical results even for rather short memory 
intervals with c < 1. In particular, the estimate for the 
fermion growth rate will be insensitive to the value of 
c. Furthermore, we Taylor expand the integrand in (|54[) 
around the latest time of the memory integrals, where 
we keep only the lowest order in the expansion. Us- 
ing the antisymmetry in time as described by (|23[) and 



(|24| the equal-time correlation functions Fy(t, t; |p|) and 
Pv{t,t; |p|) vanish, such that the corresponding terms in 
dSU) can be neglected to lowest order. Using (f3"2"| we 
note that py(t,t;\p\) = i for all momenta and times. 
Applying these approximations, the memory expanded 
evolution equation reads 



d t F v (t,t;\p\) + K(t)F v (t,t;\p\) ~ 0. 



(55) 



where we define 
K(t) 



/ F (M;|p-q| 

70 Jq 



To 



(56) 



9 2(7rt)3/2 e 



In order to perform the integral in (|56|) . we used (|53|) 
and the fact that the contributions from high momenta 
are exponentially suppressed such that it may be well 
approximated by a Gaussian integral. We note that 
all of the above analytical approximations respect chi- 
ral symmetry, which is also discussed in more detail in 
Appendix [B] Wc emphasize that the above approxima- 
tions are not valid for 27oi <SC 1 since F^t, t; \p\) is not 
large initially and one observes that K (t) is not defined 
at t = 0. However, we will see that even at early times 
the description is still reasonable and one can solve ([53]) 
with the initial condition Fv(t,t; |p|)|t = o+ = 1/2 in ac- 
cordance with (|3"5l) . 

The solution of the equation of motion can be deter- 
mined analytically and is given by 



F v (t,t;\p\) = 
-exp \-g 2 ae 2lat 



2daw(^27 £) 



V27o^ 



(57) 



with the dimcnsionlcss constant 

_ A 7oV2c 
a - ^3/2 ■ 

The Dawson function appearing in (|57[) is defined as 



(58) 



daw(x) = e 



dy e v 



(59) 



which is closely related to the error function [48j , It is 
an odd function and its derivative is daw (x) = 1 — 
2xd&w(x) such that it grows linearly at x = 0. Its max- 
imum occurs at xq ~ 0.924 and daw(xo) — 0.541. It has 
the asymptotic series 



(60) 



with which the expression (|57[) can be simplified for not 
too early times by using 



2daw(^2 7o t) - 



1 



1 



y^t 2(2 70 t)3/ 2 



(61) 
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FIG. 10: Evolution of n^/j (t, |p|) as a function of time for 
different momenta and gV2 = 10 -2 . The boson self-coupling 
is A = 10~ J . The solid line represents the analytic solution 
and the dotted lines show the results from simulations. 



FIG. 11: Evolution of n^(t, |p|) as a function of time. The 
analytical solution (solid lines) is compared to the result from 
numerical simulations (dotted lines) for different values of the 
Yukawa coupling g. For the numerical simulations the lowest 
available momentum |p] = 0.13m is chosen. 



One finds that (|6ip provides a rather good description 
already for times 270 1 > 0(1). 

The behavior of the fermion occupation number 
n$(t, |p|) is obtained by plugging the solution (|57p into 
(|4TI|l . With A 7o ~ 1/4 the constant flS5J) is approx- 
imately a ~ 0.05c. We typically find best agreement 
with the numerical solution for c ~ 0.48, which is fur- 
ther discussed below. A very characteristic feature of 
the solution ([5"7| is that it predicts a growth rate which 
is independent of momenta. In Fig. [TO] simulation re- 
sults for n^(t, |p|) for different momenta are compared 
to the corresponding analytical solution for g^/2 = 10~ 2 
and A = 10 -5 . There are deviations at early times for 
the reasons mentioned above. In particular the asymp- 
totic growth rate is accurately reproduced and indeed it 
is approximately momentum independent, which can also 
be seen from Fig. [5] The quantitative agreement of the 
analytical and full numerical solution is good for small 
momenta as expected, with sizeable corrections already 
for |p| ~ m. 

In the following we will discuss the solution (|5T)) for 
two characteristic time ranges (I) 27o£ > 0(1) and (II) 
27o< 3> 1 for weak coupling g <C 1. The r.h.s. of ([57]) 
contains two exponential functions. For the regime (I) 
the argument of the first exponential appearing in ([57)) is 
small for weak coupling. Employing a Taylor expansion 
it behaves as 



F v {t,t;\p\) ~ 



2daw(V27oi) 



V27oi 



(62) 



Accordingly, the fermion occupation number n^(t, |p|) 
given by (|40|) grows exponentially with approximately 
the same rate as the maximum primary growth rate for 



the bosons: 

n*(t, |p|) ^V^< 
g 2 a 



2daw( 1 /27o£) - 



4(27ot) 3 /2 



,2 70* 



(63) 



To arrive at the second approximate equality we used 
(|6ip . The exponential growth is parametrically slowly 
modified by the factor ~ t~ 3 / 2 , reducing somewhat the 
growth. By fitting the time dependence to just an ex- 
ponential in the respective regimes (sec e.g. Fig. 1 1 1 j) one 
may infer an approximate effective growth rate for the 
fermion modes of about 7^ w O.8570, i.e. smaller than 
the maximum primary boson growth rate. 

Because of the exponential growth the employed Taylor 
expansion that we used to obtain ([6"2"[) or ([6"3"[) becomes 
invalid for the regime (II). Without the Taylor expansion, 
using (JHTJ) in ([57| . one finds 



r 1 is 1 1 

»w> IpI) - 2 ~ 2 exp 



g 2 a 



2(2 7o t)3/2 



p2 70* 



(64) 



The asymptotic value n^(t, |p|) = 1/2 is approached 
from below with a double-exponential time dependence. 
This corresponds to a very abrupt shutoff of the exponen- 
tial growth period, as can also be seen in the behavior of 
the numerical solutions in Figs. 151 and [TOl 

The characteristic time for the end of the exponen- 
tial growth can be estimated to be the time when the 
argument of the first exponential in the solution (|57[) is 
of order one. Approximately this can be obtained from 
([M| by equating the absolute value of the argument of 
the first exponential to one, which yields 

h - ^~ ln ( 4-] + t- 1x1 (2 70 • (65) 

270 \g 2 aj 470 
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FIG. 12: Evolution of n^(t, |p|) as a function of time for 
given Yukawa coupling g and different values for the boson 
self-coupling A as indicated in the graph. The solid line cor- 
responds to the analytical solution and the dotted lines to the 
results from simulations. 



FIG. 13: Fermion occupation numbers as a function of mo- 
mentum at fixed time t — 90/m after the exponential growth 
period. Compared are simulations with different A and g but 
equal £ = g 2 /X. 



Though the growth rate does not depend on the choice of 
c appearing in the dimcnsionlcss constant a given by (|58[) . 
we note that the time depends logarithmically on it. 
Using again c ~ 0.48 we find t 4 , ~ 11.8 for g\/2 = 0.001, 
t v , ~ 9.2 for gy/2 = 0.01, and ~ 6.7 for ,gV2 = 0.1. 

In Fig.[TT]we show a comparison of the analytical result 
(|57p and the numerical simulation for different Yukawa 
couplings. The simulations are for the fermion mode with 
|p| = 0.13m (our smallest momentum) and a boson self- 
coupling of A = 10" 5 . One observes from Fig. [TTI that 
the abrupt shutoff of the exponential growth is well re- 
produced for a wide range of Yukawa couplings g, and in 
these cases our above estimates for agree rather well 
with the numerical findings. However, sizeable deviations 
occur for g\[2 < 10~ 3 with A = 10~ 5 . Correspondingly, 
Fig- HH shows results for fixed g and different values of A. 
The same statements apply as before, and similar devi- 
ations after the exponential growth terminates occur for 
A > 10 -3 with q-y/2 = 10~ 2 . One observes that in both 
cases the ratio [52| 

i = y (66) 

at which significant corrections appear in the analytical 
description of the later stages of the fermion evolution 
is about 2£ ~ 0.1. The limitation comes from the fact 
that we use the exponentially growing solution l|53p at all 
times, which is only correct before the growth of boson 
modes terminates, i.e. for times t < given by (|48[) . It 
is remarkable that according to Figs. [TT1 and [T2l this leads 
to practically no corrections to the solution (|57p also for 
later times if 2£ > 0.1. These findings can be understood 
from the fact that the characteristic time at which 
the growth of fermion occupation numbers terminates is 
earlier than tj, in these cases. We use (|B"51) together with 
(|4"5|) to estimate for which range of couplings one finds 



<0 ^t^. Without further approximations this condition 
translates into the compact expression £ > 1/ (8c) , which 
is in reasonable agreement with the numerical findings as 
explained above. 

It is a general fact that the ratio (|66|) plays an im- 
portant role for the fermion dynamics for t > 4a. This 
can be understood as follows. According to the discus- 
sion in Sec. IHI Al at the time tj, the parametric coupling- 
dependence of the boson correlation function becomes 
Fff, ~ 0(1/ A). The boson-fermion loop contribution to 
the fermion propagators shown in Fig. [9] is proportional 
to g 2 . Since for the dominant term the boson propaga- 
tor line appearing in the loop is associated to F^, this 
leads to the parametric g 2 / X or £ dependence. In order 
to verify this we have plotted in Fig. [13] fermion occupa- 
tion number distributions for different values of g and A. 
Compared are different simulations with same ratio £ for 
2£ = 0.1, 0.01 and 0.001 at fixed time t = 90/m after the 
exponential growth period. The agreement is remark- 
ably good. We emphasize that at early times they differ 
in general for the same £. In particular, the time itself 
depends on A and is insensitive to g for the considered 
weak couplings (see Sec. IIII A"| . 



B. Emergence of power-law behavior 

We have seen in Sec. IIII Bl that boson occupation num- 
bers approach a nonthermal scaling solution ~ p _K iR 
with occupation number exponent Kir, ~ 4 for low- 
momentum modes after the exponential growth period. 
The corresponding nonthermal infrared fixed point has 
the dramatic consequence that a diverging time scale 
exists far from equilibrium, which can prevent or sub- 
stantially delay thermalization due to critical slowing 
down (l3l , [ItJ • Because fermion field degrees of freedom 
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obey the Pauli exclusion principle their occupation num- 
bers cannot become large (53[. This preempts infrared 
scaling solutions, which diverge in the infrared. As a 
consequence, the phenomenon of critical slowing down 
is expected to be absent for fermions. Low-momentum 
fermion modes may thus approach a thermal distribution 
much faster than bosons in this case, which we observed 
already in Sec. IIV Al and discuss in more detail below. 

In contrast to infrared scaling solutions, power-law be- 
havior for higher momentum modes, where occupation 
numbers are not large, is of course not excluded for 
fermions. In the right column of Fig. [5J which shows 
the evolution of n^(t, |p|), the fermion occupation num- 
bers exhibit a power-law regime with exponent two for 
momenta |p| > m. We have verified that the very sharp 
fall-off at later times of the distribution close to the cut- 
off, i.e. |p| w 10m in Fig. is removed if we enlarge 
the momentum cutoff and the power-law behavior ex- 
tends to larger momenta accordingly. One observes from 
Fig. [5] that the emergence of power-law behavior happens 
quite early, i.e. around the time t m 3/m and it is well es- 
tablished already around t ~ 5/m. This points out that 
the origin of the fermion power-law is distinct from the 
scaling behavior of the bosons discussed in Sec. IIII Bl 

In the following we will explain the fermion power-law 
regime. From Fig. [51 in particular from the upper graph 
of that figure, one observes that the growth of higher 
momentum modes sets in later than for lower momenta. 
Since the growth rate is approximately momentum in- 
dependent (see Sec. IIV A[) the retardation in the initial 
growth time translates into a momentum dependent am- 
plitude 



A(\P\) 



n^(t, |p|) 



O 2 7o t 



(67) 



The residual time-dependence of this amplitude is para- 
metrically slow compared to the dominant exponential 
growth behavior of n^{t, |p|). The observed power-law 
with exponent two then corresponds to A(|p|) ~ |p| -2 
in the respective momentum range with |p| > to. A fur- 
ther important property, which may be observed from 
Fig. [8J but better from the right column of Fig. [5l is 
that the time for the end of the exponential growth pe- 
riod becomes approximately independent of momentum 
for not too small momenta. As a consequence, part of 
the momentum-dependence emerging at early times be- 
comes "frozen- in" in time during the subsequent quasi- 
stationary evolution. In order to quantify these state- 
ments, we repeat the analytical calculation of Sec. IIV Al 
with the main difference that we take retardation or 
memory effects into account. To keep the discussion an- 
alytically tractable, our approximations will describe the 
dominant exponential time-dependence while subleading 
powers in time are discarded in contrast to what was done 
in Sec. IIV Al This does not affect the relevant momen- 
tum dependence of -A(|p|) but overestimates the overall 
size of the amplitude. 

Starting point of the analysis is again the equation 



of motion for Fy(t, t; |p|) presented in ([54]) . Approxi- 
mating the boson two-point function F^it, t'\ |p|) by its 
primary-growth behavior (|44[) and applying a memory 
expansion as in Sec. IIV Al would lead to the local differen- 
tial equation with respect to momentum (]55p . such that 
the momentum-independent solution (|57p is obtained for 
the employed initial conditions. While this is a good 
approximation for the considered low-momentum behav- 
ior with |p| < to in Sec. IIV Al here we concentrate on 
larger momenta. In order to make analytical progress, 
a perturbative estimate of the fermion correlators under 
the integral (p>4"]) is employed for small g. This amounts 
to replacing the respective correlators by their free-field 
solutions: 



1 



cos(|p|(t - t 1 ' 



- l -sm(\p\{t 
sin(|p|(t-0 
icos(|p|(i-i" 



-t" 



(68) 
(69) 
(70) 
(71) 



F v (t,t";\p\) 

F v (tX;\p\) 
pv{t, t";\p\) 
P°v(tX;\p\) 



With these expressions the equation of motion (j54|) be- 
comes 

d t F v (t,t;\p\) ~ 

- 2.g 2 | t dt"J F^t"; IP - q|) cos(2|p|(t - i")) . 

(72) 

Using accordingly the lowest-order behavior (j44| for the 
scalar correlator, where F$(t,t";\p\ > to) vanishes or 
may be approximated by a Gaussian as in (|56[) . the mo- 
mentum integral in (|72p can be directly evaluated. How- 
ever, the following time-integration would lead to expres- 
sions which are rather inconvenient to handle. Therefore, 
we may simplify the approximation by using 7(p) ~ 70 
in the momentum integral in (]72p . such that 



F^(t,t" 



6tt 2 



(73) 



Though this captures the dominant exponential time- 
dependence, comparison to (|56p shows that this ne- 
glects a parametrically slow (t + £")~ 3 / 2 -correction. Wc 
have verified that this does not affect our main conclu- 
sions about momentum-dependencies, however, it over- 
estimates the overall growth. Pursuing with (|73p in the 
equation of motion (|72p the memory integral can be eval- 
uated and gives the compact expression 



d t F v (t,t;\p\ 



7 e 27t>t -e 70 * 



3tt 2 

7 cos(2|p|i)-2|p|sin(2|p|i) 



4p 



7o 2 



(74) 
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FIG. 14: Comparison between Fermi-Dirac distribution (thin) 
and simulation of the fermion occupation number at t = 90/m 
(thick) for 2£ = 1 (solid line) and 2£ = 0.01 (dotted line). 



Using the initial condition (|55|) this can be straightfor- 
wardly integrated and one obtains for the fermion occu- 
pation number (|40p : 



n^(t, |p|) 



g 2 A o7 ge 2 ^*-2e^ t cos(2|p|t) 

6tt 2 4|p7 
«? 2 A o7o 3 e 2 ^ 4 



7o 2 



6tt 2 4|p| : 



7o 2 



(75) 



While the first line in (|75| represents a solution of (|74|) 
without further approximations, for the second line we 
exploit the fact that the dominant exponential term 
quickly takes over. Thus the oscillatory behavior be- 
comes suppressed. We conclude from (|75|) that the am- 
plitude 



A(\p\) 



1 



4|p| 2 +7o 2 



(76) 



which is in reasonable agreement with the behavior of 
the numerical solution displayed in the right column of 
Fig. [5] In particular, with 70 = m it describes the ob- 
served power-law exponent two for |p| large enough com- 
pared to to. We emphasize that the above perturbative 
evaluation using ([6"8"]) - ([7"Tj) cannot describe the termina- 
tion of the exponential growth period, in contrast to the 
self-consistent estimate of Sec. HVAl leading to (|5"T|) . 

We have mentioned above that low-momentum fermion 
modes may approach a thermal distribution much faster 
than bosons. More precisely, the observed quasi- 
stationary solution then shares the property of the equi- 
librium solution that n^,(|p| = 0) = 1/2 in the infrared 
[5H ]. Following the discussion of Sec. IIV Al this is ex- 
pected for couplings such that the ratio (|66|) is 2£ > 0.1. 
In this case the fermion occupation number growth ter- 
minates when the thermal equilibrium distribution is ap- 
proached in the infrared. In order to quantify this aspect 



we compare in Fig. [14] a Fermi-Dirac distribution with 
the occupation number n^,(t, |p|) in the quasi-stationary 
regime at t = 90/m. The results for 2£ = 1 are obtained 
from gV2 = 0.01 and A = 10~ 4 , whereas 2£ = 0.01 re- 
sults from g\[2 = 0.01 and A = 0.01. More precisely, 
we compare to the thermal distribution of free massless 
particles with zero net charge, 



"fd(IpI) = 



1 



elpl/ T + l 



(77) 



where the temperature is estimated by equating the re- 
spective energies obtained in a quasi-particle picture, 



d 3 p 
(2^)3 



|p|n FD (|p|) 



d 3 p 
(2^)3 



|p|Mt,|p|). (78) 



Fig. [14] indeed shows that for 2£ = 1 the quasi- 
stationary occupation number follows the thermal dis- 
tribution rather closely up to momenta |p| ~ to. This 
has to be compared to the behavior of the bosons, which 
still evolve in the infrared orders of magnitude in time 
longer towards the nonthermal fixed point distribution 
ritj, (t, |p|) ~ I p> I 4 , which is far from thermal equilib- 
rium (see Sec. UTlB|) . Fig. [Til shows that for 2£ = 0.01 
the exponential fermion occupation number growth stops 
before tifdGpI = 0) = 1/2 is reached. Accordingly, 
the deviations from the thermal distribution are larger 
in the infrared and also the differences at higher mo- 
menta due to the nonthermal power-law behavior are 
even more pronounced. Phenomenologically a major as- 
pect is that the power-law leads to a substantial excess 
of high-momentum particles compared to thermally pro- 
duced particles. 



V. CONCLUSIONS 

We have studied noncquilibrium fermion production 
in a Yukawa-type model in 3 + 1 dimensions following 
a spinodal/tachyonic instability. The approximation is 
based on a 1/N expansion to NLO of the 2PI effective 
action out of equilibrium, including the two-loop boson- 
fermion fluctuations displayed in Fig. [JJ This resums 
an infinite series of scattering processes and takes into 
account memory effects. 

These quantum corrections turned out to be crucial for 
the phenomenon of instability-induced fermion produc- 
tion observed in this work. This production mechanism 
is based on the exponential growth of long-wavelength 
boson occupation numbers following a noncquilibrium 
instability. A major characteristic property is that the 
amplification of fermion occupation numbers is approx- 
imately momentum-independent and proceeds with the 
maximum primary growth rate of boson modes. 

Depending on the ratio £ of the square of the Yukawa 
coupling and the boson self-coupling, the end of the ex- 
ponential amplification of fermion modes is either deter- 
mined by the end of the boson growth period (2£ < 0.1), 
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or by reaching the Fermi-Dirac distribution in the in- 
frared for momenta |p| < m (2£ > 0.1). The latter leads 
to very fast thcrmalization of low- momentum modes. 
This finding has to be confronted with the extremely slow 
thermalization of low- momentum bosons, which are still 
far from equilibrium at this time. 

The slowing-down of the boson dynamics after the ex- 
ponential growth period is shown to be associated to 
the approach to a nonthermal infrared fixed point. It 
is characterized by scaling of occupation number modes 
"^(IpI) ~ | P> | 4 , which lead to strongly enhanced fluc- 
tuations in the infrared as compared to a thermal low- 
momentum distribution ~ |p| _1 - It was pointed out in 
Ref. [l3j that such a fixed point is approached in scalar 
models after a parametric resonance instability, which 
can prevent or substantially delay thermalization. Here 
we have demonstrated that a spinodal/tachyonic insta- 
bility leads to scaling with the same occupation number 
exponent, which suggests that this is a universal phe- 
nomenon of nonequilibrium instability dynamics. We 
note that for the spinodal/tachyonic instability and the 
considered range of parameters we seem not to observe 
the expected phenomenon of perturbative Kolmogorov 
wave turbulence at higher momenta Instead boson 
modes with momenta |p| > m exhibit thermal properties 
at this stage. 

Fermion occupation number modes are bounded from 
above, which preempts infrared scaling solutions. We 
have explicitly verified that the absence of a scaling 
regime for fermions at low momenta does not affect the 
infrared scaling behavior observed for bosons. In gen- 
eral we find that the boson dynamics is rather insensi- 
tive to the presence of fermions for not too large Yukawa 
coupling g < 1, which is expected because of the com- 
parably large boson occupation numbers. Despite the 
absence of an infrared scaling regime for fermions, we 
find a power-law behavior with exponent two for higher 
momenta |p| > m. Its origin is very distinct from the 
physics underlying nonthermal fixed points. The latter 
appear at late times and can be described as stationary 
points of time evolution equations that are independent 
of the initial conditions. In contrast, we could show that 
the observed fermion power-law behavior emerges at very 
early times and is a consequence of a retardation or mem- 
ory effect. This power- law leads to an excess of highly 
energetic particles compared to thermal distributions. 

It is remarkable that for the considered model these 
nonlinear out-of-equilibrium phenomena can be studied 
numerically - and to a large extent analytically - using 
2PI effective action techniques directly in quantum field 
theory. An important improvement would be to include 
a macroscopic field into our analysis, in order to estimate 
the importance of quantum corrections to frequently 
employed approximations based on the Dirac equation 
in the presence of classical fields. Understanding the 
impact of quantum fluctuations is an important step 
towards more realistic phcnomcnological applications in 
early universe cosmology as well as QCD in the context 
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FIG. 15: Time evolution of the unequal-time correlation func- 
tion p$(t, £'; |p|) for different t' with couplings and momentum 
as indicated in the figure. For each t' the analytical expression 
(|A1|) is plotted for comparison. 



of collision experiments of heavy nuclei. 
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APPENDIX A: TIME EVOLUTION OF THE 
SPECTRAL FUNCTION 

In this Appendix we consider the evolution of the 
unequal-time boson spectral function p<j,(t, t'; |p|) and 
statistical function F,p(t,t'; |p|). We will use this to ex- 
plain in more detail some estimates in the main text 
based on the dominance of contributions coming from the 
equal-time correlator F$(t,t; |p|). According to the 
equal-time spectral function p^(t,t; |p|) vanishes identi- 
cally. The time evolution of F^,(t, t'; |p|) in the linear 
regime can be approximated by (|44[) . Similar considera- 
tions for the spectral function lead to 

^(t, 4 ';|p|) = i( e 7(P)(*-t')_ e -7(p)(t-*')) , (Al) 

with 7(p) defined in In Fig. [T5] we plot results 

from simulations together with the corresponding ana- 
lytical expression (|A1[) for couplings and momentum as 
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FIG. 16: Comparison of the time evolution of the statistical 
and spectral function together with the analytical expressions 
discussed in the text. Couplings and the chosen momentum 
are indicated in the figure. 



indicated in the figure. The agreement before the ex- 
ponential amplification terminates is remarkable for low- 
momentum modes. For momenta |p| > m the approxi- 
mation becomes worse. With (|AT]) and flU]) for the sta- 
tistical function one observes that for not too early times 
F^tft; |p|) dominates over F^(f, i'; |p|) and p^t^t'; |p|). 
Moreover one can find in this time regime the relation 



1 



i^(t,0;|p|)~ -^(t,0;|p 



(A2) 



which holds also for large times. In Fig. [Tj)] the corre- 
lation functions F<f,(t,t; |p|), F^i, 0; |p|) and p^(i,0; |p|) 
are plotted together with and the dominant part of 
(|A1[) with a factor 1/2. This confirms (|A2[) and the dom- 
inance of F,f,(t,t; |p|) is well established for not too early 
times. We note that for t' > the correlator F^(i, t'; |p|) 
becomes larger than F$(t, 0; |p|), whereas p$(t, t'\ |p|) be- 
comes smaller than p<t,(t, 0; |p|) as Fig. [T51 shows. These 
findings justify the neglection of terms proportional to 
p<l>{t,t';\p\) in the self-energy contributions in (jTTJ) , as 
was done for the analytical estimates in Sec. IIVI 



APPENDIX B: DECOMPOSITION OF FERMION 
TWO-POINT FUNCTIONS 

In this Appendix we present the expressions of the 
statistical- and spectral-function of the fermion two-point 
function in the chiral basis. The statistical function is de- 
fined by (we omit Dirac indices) 



(Bl) 



where the expectation value has to be taken with respect 
to the initial density matrix and = ip^x)^. In the 



chiral basis we can write ip(x) = i/jl(x) + ipn(x). Where 



^l(x) 



Xl(x) 




and ipn(x) 





Xr(x) 



(B2) 



The Xl/r(x) are two-component Weyl spinors. This de- 
composition is done by projections (Pl and Pr) given 

as 

if>(x) = P L ^{x) + Pr^(x) 

ee 1(1- 75 )tP( x ) + 1(1 +75 )iP( x ). (B3) 

The 7-matrices in this basis are given by 



7o 

75 



l 2 x2 
1 2 X2 

-l2xa 

l 2 x2 



CTi 

-d 



(B4) 



Using this notation we get an expression for the statis- 
tical function in terms of left and right handed spinors. 
Here we merely give the result of the vector components 
relevant in the analysis given in the main text. The ex- 
pression for the temporal component is 



1 



tr 



+tr 



Xl(x),Xl(v) 



(B5) 



and for the spatial components one finds 
1 



*V(x,y) 



-tr er 



8 



tr 



<r( {xr{x),Xr(v)} 



[xL(x),xl(y)} 

The spectral function is defined by 

p(x,y) = i( {ip(x),ip{y)} 



(B6) 



(B7) 



Similarly as above, the temporal component of the spec- 
tral function is 



tr 



{xr(x),X r (v)] 



+tr \({xL(x),xi(y)} 

and the spatial components are given by 



(B8) 



tr a 



-tr(cr( XL(x),x{{y) 



XR.(x),x R (y) 



(B9) 
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APPENDIX C: NUMERICAL METHOD 

Our numerical calculations for the equations of motion 
derived from the 2PI effective action are carried out on a 
single personal computer with 8Gb RAM. The discretiza- 
tion is performed on the level of the equations of motion. 
We discretize two time dimensions and use spherical sym- 
metry for the three-dimensional momentum integrals, so 
that only the absolute value of the discretized momenta 
needs to remain in memory. The angular part of the in- 
tegrals can be done analytically in each case because the 
arising convolutions in momentum space are given as a 
product of two functions in coordinate space. The tran- 
sition between coordinate and momentum space is imple- 
mented using the FFTW library ( |http://www.fftw.org/p 
for one-dimensional Fourier transforms. For the time di- 
mension we use up to 1850 time steps with a width of 
dt < 0.05, depending on the coupling constants. The 
momentum grid has N = 80 discretization points with 
an equidistant spacing of dp = tt/(N a s ), where a s = 0.3 
is the spatial grid spacing. We use a leap-frog- typ e al- 
gorithm following earlier work, for details see [iji l35j |. 
With this algorithm the fermion sector has a discretiza- 



tion in time with a time-step width of twice the one in 
the boson sector. All of the derivatives and integrals are 
calculated with the Eulcr method and the trapezoid rule, 
respectively. 

In the case of the classical-statistical simulations for 
bosonic fields (see Figs. |6]and[2|), we use the standard 
hypcrcubic lattice in 3+1 dimensions. The employed 
Laplace operator in the classical field equation of motion 
is of fourth order and combined with a finite-difference 
second-order time derivative. The field initial conditions 
are generated to match the Gaussian initial correlators 
for the quantum 2PI NLO calculation and a 'quench' 
is implemented by changing the sign of the boson mass 
within the two time-steps necessary to set up a second or- 
der time evolution. Averages for two-point functions are 
obtained by summing over all possible redundant com- 
binations of distances or momenta, respectively. The 
necessary transition between coordinate and momentum 
space is implemented using the FFTW library for three- 
dimensional Fourier transforms. Our maximum possible 
lattice size on a 8Gb memory machine is 460 points along 
each spatial axis with the corresponding spacing a s = 0.3 
whereas the each temporal step corresponds to dt = 0.05. 
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